home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_10_08 / 1008038a < prev    next >
Encoding:
Text File  |  1992-05-28  |  16.4 KB  |  422 lines
  1. /* RAY_RAD - A Very Simple Ray-Traced Radiosity Renderer        */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <math.h>
  6.  
  7. #define FLOOR           0       /* X-Y plane, Z = 0.0           */
  8. #define SOUTH_WALL      1       /* X-Z plane, Y = 0.0           */
  9. #define EAST_WALL       2       /* Y-Z plane, X = 8.0           */
  10. #define NORTH_WALL      3       /* X-Z plane, Y = 8.0           */
  11. #define WEST_WALL       4       /* X-Z plane, X = 0.0           */
  12. #define MAX_RAYS        25000   /* Maximum number of rays       */
  13. #define MAX_VAL         1.0e+6  /* Maximum floating point value */
  14. #define MIN_VAL         1.0e-6  /* Minimum floating point value */
  15. #define NUM_LOOP        20      /* Number of iterations         */
  16. #define NUM_SURF        10      /* Number of surfaces           */
  17. #define PI              3.141592654
  18. #define DRAND()         ((double) rand() / (double) RAND_MAX)
  19.  
  20. typedef struct          /* 3-D point co-ordinates               */
  21. {
  22.   double x;             /* X-axis co-ordinate                   */
  23.   double y;             /* Y-axis co-ordinate                   */
  24.   double z;             /* Z-axis co-ordinate                   */
  25. } PT_3D;
  26.  
  27. typedef struct          /* 3-D ray                              */
  28. {
  29.   PT_3D org;            /* Origin                               */
  30.   PT_3D dir;            /* Direction                            */
  31. } RAY;
  32.  
  33. typedef struct          /* Element                              */
  34. {
  35.   double total;         /* Total flux                           */
  36.   double unsent;        /* Unsent flux                          */
  37. } ELEM;
  38.  
  39. typedef struct          /* Surface                              */
  40. {
  41.   PT_3D vertex[4];      /* Vertex array                         */
  42.   int num_row;          /* Number of element rows               */
  43.   int num_col;          /* Number of element columns            */
  44.   double col_dim;       /* Element column dimension             */
  45.   double row_dim;       /* Element row dimension                */
  46.   double emittance;     /* Emittance                            */
  47.   double reflectance;   /* Reflectance                          */
  48.   ELEM *elemp;          /* Element array pointer                */
  49. } SURF;
  50.  
  51. static int send_flux(int, int, int);
  52. static void display_surface(void);
  53. static void find_element(RAY *, int *, int *, int *);
  54. static void global_to_local(int, PT_3D *, PT_3D *);
  55. static void local_to_global(int, RAY *, RAY *);
  56. static void select_ray(int, int, int, RAY *);
  57.  
  58. static double MaxEmittance = 0.0;       /* Maximum emittance    */
  59. static SURF RoomSurf[NUM_SURF] = {      /* Room surfaces        */
  60.   { { { 0.0, 0.0, 0.0 }, { 8.0, 0.0, 0.0 },     /* Floor        */
  61.     { 8.0, 8.0, 0.0 }, { 0.0, 8.0, 0.0 } }, 8, 8, 1.0, 1.0,
  62.        0.0, 0.2, NULL },
  63.   { { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 8.0 },     /* South wall   */
  64.     { 8.0, 0.0, 8.0 }, { 8.0, 0.0, 0.0 } }, 8, 8, 1.0, 1.0,
  65.        0.0, 0.5, NULL },
  66.   { { { 8.0, 0.0, 0.0 }, { 8.0, 0.0, 8.0 },     /* East wall    */
  67.     { 8.0, 8.0, 8.0 }, { 8.0, 8.0, 0.0 } }, 8, 8, 1.0, 1.0,
  68.        0.0, 0.5, NULL },
  69.   { { { 8.0, 8.0, 0.0 }, { 8.0, 8.0, 8.0 },     /* North wall   */
  70.     { 0.0, 8.0, 8.0 }, { 0.0, 8.0, 0.0 } }, 8, 8, 1.0, 1.0,
  71.        0.0, 0.5, NULL },
  72.   { { { 0.0, 8.0, 0.0 }, { 0.0, 8.0, 8.0 },     /* West wall    */
  73.     { 0.0, 0.0, 8.0 }, { 0.0, 0.0, 0.0 } }, 8, 8, 1.0, 1.0,
  74.        0.0, 0.5, NULL },
  75.   { { { 0.0, 0.0, 8.0 }, { 0.0, 3.0, 8.0 },     /* South ceil.  */
  76.     { 8.0, 3.0, 8.0 }, { 8.0, 0.0, 8.0 } }, 8, 3, 1.0, 1.0,
  77.        0.0, 0.8, NULL },
  78.   { { { 5.0, 3.0, 8.0 }, { 5.0, 5.0, 8.0 },     /* East ceiling */
  79.     { 8.0, 5.0, 8.0 }, { 8.0, 3.0, 8.0 } }, 3, 2, 1.0, 1.0,
  80.        0.0, 0.8, NULL },
  81.   { { { 0.0, 5.0, 8.0 }, { 0.0, 8.0, 8.0 },     /* North ceil.  */
  82.     { 8.0, 8.0, 8.0 }, { 8.0, 5.0, 8.0 } }, 8, 3, 1.0, 1.0,
  83.        0.0, 0.8, NULL },
  84.   { { { 0.0, 3.0, 8.0 }, { 0.0, 5.0, 8.0 },     /* West ceiling */
  85.     { 3.0, 5.0, 8.0 }, { 3.0, 3.0, 8.0 } }, 3, 2, 1.0, 1.0,
  86.        0.0, 0.8, NULL },
  87.   { { { 3.0, 3.0, 8.0 }, { 3.0, 5.0, 8.0 },     /* Ceil. light  */
  88.     { 5.0, 5.0, 8.0 }, { 5.0, 3.0, 8.0 } }, 2, 2, 1.0, 1.0,
  89.     5000.0, 0.8, NULL } };
  90.  
  91. int main(void)
  92. {
  93.   int col;              /* Column counter                       */
  94.   int elem;             /* Element counter                      */
  95.   int maximum;          /* Maximum number of rays               */
  96.   int num_elem;         /* Number of elements                   */
  97.   int num_rays;         /* Number of rays                       */
  98.   int row;              /* Row counter                          */
  99.   int surf;             /* Surface counter                      */
  100.   ELEM *elemp;          /* Element pointer                      */
  101.   SURF *surfp;          /* Surface pointer                      */
  102.  
  103.   for (surf = 0; surf < NUM_SURF; surf++) {
  104.     /* Instantiate elements                                     */
  105.     surfp = &(RoomSurf[surf]);
  106.     num_elem = surfp->num_row * surfp->num_col;
  107.     if ((surfp->elemp = (ELEM *) malloc(num_elem * sizeof(ELEM)))
  108.         == NULL) {
  109.       fputs("Out of memory!\n", stderr);
  110.       return (2);
  111.     }
  112.  
  113.     elemp = surfp->elemp;
  114.     for (elem = 0; elem < num_elem; elem++, elemp++)
  115.       elemp->total = elemp->unsent = surfp->emittance;
  116.  
  117.     if (surfp->emittance > MaxEmittance)
  118.       MaxEmittance = surfp->emittance;
  119.   }
  120.  
  121.   do {
  122.     /* Distribute flux between elements                         */
  123.     num_rays = 0;
  124.     for (surf = 0; surf < NUM_SURF; surf++) {
  125.       surfp = &(RoomSurf[surf]);
  126.       for (row = 0; row < surfp->num_row; row++)
  127.         for (col = 0; col < surfp->num_col; col++) {
  128.           if ((maximum = send_flux(surf, row, col)) > num_rays)
  129.             num_rays = maximum;
  130.         }
  131.     }
  132.     display_surface();          /* Display intermediate results */
  133.   } while (num_rays > 0);       /* Repeat until convergence     */
  134.  
  135.   for (surf = 0; surf < NUM_SURF; surf++)       /* Free memory  */
  136.     free(RoomSurf[surf].elemp);
  137.  
  138.   return (0);
  139. }
  140.  
  141. static int send_flux (  /* Send flux to other elements          */
  142.   int s_surf,           /* Sending surface index                */
  143.   int s_row,            /* Sending element row                  */
  144.   int s_col )           /* Sending element column               */
  145. {
  146.   int h_col;            /* Hit element column                   */
  147.   int h_row;            /* Hit element row                      */
  148.   int h_surf;           /* Hit surface index                    */
  149.   int num_rays;         /* Number of rays                       */
  150.   int ray;              /* Ray counter                          */
  151.   double inc_flux;      /* Incident (ray) flux                  */
  152.   double ref_flux;      /* Reflected flux                       */
  153.   RAY shoot;            /* Shooting ray                         */
  154.   ELEM *h_elemp;        /* Hit element pointer                  */
  155.   ELEM *s_elemp;        /* Sending element pointer              */
  156.   SURF *h_surfp;        /* Hit surface pointer                  */
  157.   SURF *s_surfp;        /* Sending surface pointer              */
  158.  
  159.   s_surfp = &(RoomSurf[s_surf]);
  160.   s_elemp = s_surfp->elemp + s_row * s_surfp->num_col + s_col;
  161.  
  162.   /* Number of rays to shoot proportional to unsent flux        */
  163.   if ((num_rays = (int) (MAX_RAYS * s_elemp->unsent /
  164.       MaxEmittance)) == 0)
  165.     return (0);
  166.  
  167.   inc_flux = s_elemp->unsent / num_rays;    /* Get ray flux     */
  168.  
  169.   for (ray = 0; ray < num_rays; ray++) {    /* Distribute flux  */
  170.     select_ray(s_surf, s_row, s_col, &shoot);   /* Select ray   */
  171.  
  172.     /* Find hit surface and element                             */
  173.     find_element(&shoot, &h_surf, &h_row, &h_col);
  174.  
  175.     /* Calculate flux reflected from hit element                */
  176.     h_surfp = &(RoomSurf[h_surf]);
  177.     h_elemp = h_surfp->elemp + h_row * h_surfp->num_col + h_col;
  178.     ref_flux = h_surfp->reflectance * inc_flux;
  179.  
  180.     h_elemp->total += ref_flux;         /* Update total flux    */
  181.     h_elemp->unsent += ref_flux;        /* Update unsent flux   */
  182.   }
  183.  
  184.   s_elemp->unsent = 0.0;    /* Reset sender's unsent flux       */
  185.  
  186.   return (num_rays);
  187. }
  188.  
  189. static void select_ray (        /* Randomly select ray          */
  190.   int surf,     /* Surface index                                */
  191.   int row,      /* Element row                                  */
  192.   int col,      /* Element column                               */
  193.   RAY *rayp )   /* Ray pointer                                  */
  194. {
  195.   double horz;          /* Ray direction local horizontal angle */
  196.   double vert;          /* Ray direction local vertical angle   */
  197.   RAY local;            /* Ray local co-ordinates               */
  198.   SURF *surfp;          /* Surface pointer                      */
  199.  
  200.   surfp = &(RoomSurf[surf]);
  201.  
  202.   /* Randomly select local (surface) ray origin co-ordinates    */
  203.   /* within sending element boundary                            */
  204.   local.org.x = ((double) col + DRAND()) / surfp->col_dim;
  205.   local.org.y = ((double) row + DRAND()) / surfp->row_dim;
  206.   local.org.z = 0.0;
  207.  
  208.   /* Randomly select local (surface) ray direction co-ordinates */
  209.   /* with cosine probability distribution in vertical plane     */
  210.   horz = 2.0 * PI * DRAND();            /* Horizontal angle     */
  211.   vert = acos(sqrt(1.0 - DRAND()));     /* Vertical angle       */
  212.  
  213.   /* Convert from spherical to rectilinear co-ordinates         */
  214.   local.dir.x = sin(vert) * cos(horz);
  215.   local.dir.y = sin(vert) * sin(horz);
  216.   local.dir.z = cos(vert);
  217.  
  218.   /* Convert ray co-ordinates from local to global              */
  219.   local_to_global(surf, &local, rayp);
  220. }
  221.  
  222. static void find_element (      /* Find intersected element     */
  223.   RAY *rayp,            /* Ray pointer                          */
  224.   int *h_surfp,         /* Hit surface index pointer            */
  225.   int *h_rowp,          /* Hit element row pointer              */
  226.   int *h_colp )         /* Hit element column pointer           */
  227. {
  228.   int surf;             /* Surface counter                      */
  229.   double s = MAX_VAL;   /* Smallest parametric value            */
  230.   double t;             /* Current parametric value             */
  231.   PT_3D hit;            /* Hit point co-ordinates               */
  232.   PT_3D temp;           /* Temporary point co-ordinates         */
  233.   SURF *surfp;          /* Surface pointer                      */
  234.  
  235.   for (surf = 0; surf < NUM_SURF; surf++) {
  236.     surfp = &(RoomSurf[surf]);
  237.  
  238.     /* Determine if and where ray intersects surface plane      */
  239.     switch (surf) {
  240.       case EAST_WALL:
  241.       case WEST_WALL:
  242.         if (fabs(rayp->dir.x) > MIN_VAL) {
  243.           t = (surfp->vertex[0].x - rayp->org.x) / rayp->dir.x;
  244.           if (t > MIN_VAL && t < s) {
  245.             s = t;
  246.             *h_surfp = surf;
  247.           }
  248.         }
  249.         break;
  250.       case SOUTH_WALL:
  251.       case NORTH_WALL:
  252.         if (fabs(rayp->dir.y) > MIN_VAL) {
  253.           t = (surfp->vertex[0].y - rayp->org.y) / rayp->dir.y;
  254.           if (t > MIN_VAL && t < s) {
  255.             s = t;
  256.             *h_surfp = surf;
  257.           }
  258.         }
  259.         break;
  260.       case FLOOR:
  261.         if (fabs(rayp->dir.z) > MIN_VAL) {
  262.           t = (surfp->vertex[0].z - rayp->org.z) / rayp->dir.z;
  263.           if (t > MIN_VAL && t < s) {
  264.             s = t;
  265.             *h_surfp = surf;
  266.           }
  267.         }
  268.         break;
  269.       default:
  270.         if (fabs(rayp->dir.z) > MIN_VAL) {
  271.           t = (surfp->vertex[0].z - rayp->org.z) / rayp->dir.z;
  272.           if (t > MIN_VAL && t < s) {
  273.             /* Check for intersection inside surface boundary   */
  274.             temp.x = rayp->org.x + t * rayp->dir.x;
  275.             temp.y = rayp->org.y + t * rayp->dir.y;
  276.             if ((temp.x >= surfp->vertex[0].x && temp.x <
  277.                 surfp->vertex[3].x) && (temp.y >=
  278.                 surfp->vertex[0].y && temp.y <
  279.                 surfp->vertex[1].y)) {
  280.               s = t;
  281.               *h_surfp = surf;
  282.             }
  283.           }
  284.         }
  285.         break;
  286.     }
  287.   }
  288.  
  289.   /* Calculate hit point global co-ordinates                    */
  290.   temp.x = rayp->org.x + s * rayp->dir.x;
  291.   temp.y = rayp->org.y + s * rayp->dir.y;
  292.   temp.z = rayp->org.z + s * rayp->dir.z;
  293.  
  294.   /* Convert global to local (surface) co-ordinates             */
  295.   global_to_local(*h_surfp, &temp, &hit);
  296.  
  297.   /* Calculate hit element row and column indices               */
  298.   *h_rowp = (int) (hit.y / surfp->row_dim);
  299.   *h_colp = (int) (hit.x / surfp->col_dim);
  300. }
  301.  
  302. static void global_to_local (   /* Convert global to local      */
  303.   int surf,             /* Surface index                        */
  304.   PT_3D *globalp,       /* Global co-ordinates pointer          */
  305.   PT_3D *localp )       /* Local co-ordinates pointer           */
  306. {
  307.   PT_3D *vp;    /* Vertex pointer                               */
  308.  
  309.   vp = &(RoomSurf[surf].vertex[0]);
  310.   switch (surf) {
  311.     case FLOOR:
  312.       localp->x = globalp->x - vp->x;
  313.       localp->y = globalp->y - vp->y;
  314.       break;
  315.     case SOUTH_WALL:
  316.       localp->x = globalp->z - vp->z;
  317.       localp->y = globalp->x - vp->x;
  318.       break;
  319.     case EAST_WALL:
  320.       localp->x = globalp->z - vp->z;
  321.       localp->y = globalp->y - vp->y;
  322.       break;
  323.     case NORTH_WALL:
  324.       localp->x = globalp->z - vp->z;
  325.       localp->y = vp->x - globalp->x;
  326.       break;
  327.     case WEST_WALL:
  328.       localp->x = globalp->z - vp->z;
  329.       localp->y = vp->y - globalp->y;
  330.       break;
  331.     default:
  332.       localp->x = globalp->y - vp->y;
  333.       localp->y = globalp->x - vp->x;
  334.       break;
  335.   }
  336. }
  337.  
  338. static void local_to_global (   /* Convert ray co-ordinates     */
  339.   int surf,             /* Surface index                        */
  340.   RAY *localp,          /* Local co-ordinates pointer           */
  341.   RAY *globalp )        /* Global co-ordinates pointer          */
  342. {
  343.   PT_3D *vp;    /* Vertex pointer                               */
  344.  
  345.   vp = &(RoomSurf[surf].vertex[0]);
  346.   switch (surf) {
  347.     case FLOOR:
  348.       globalp->org.x = vp->x + localp->org.x;
  349.       globalp->org.y = vp->y + localp->org.y;
  350.       globalp->org.z = vp->z + localp->org.z;
  351.       globalp->dir.x = localp->dir.x;
  352.       globalp->dir.y = localp->dir.y;
  353.       globalp->dir.z = localp->dir.z;
  354.       break;
  355.     case SOUTH_WALL:
  356.       globalp->org.x = vp->x + localp->org.y;
  357.       globalp->org.y = vp->y + localp->org.z;
  358.       globalp->org.z = vp->z + localp->org.x;
  359.       globalp->dir.x = localp->dir.y;
  360.       globalp->dir.y = localp->dir.z;
  361.       globalp->dir.z = localp->dir.x;
  362.       break;
  363.     case EAST_WALL:
  364.       globalp->org.x = vp->x - localp->org.z;
  365.       globalp->org.y = vp->y + localp->org.y;
  366.       globalp->org.z = vp->z + localp->org.x;
  367.       globalp->dir.x = -localp->dir.z;
  368.       globalp->dir.y = localp->dir.y;
  369.       globalp->dir.z = localp->dir.x;
  370.       break;
  371.     case NORTH_WALL:
  372.       globalp->org.x = vp->x - localp->org.y;
  373.       globalp->org.y = vp->y - localp->org.z;
  374.       globalp->org.z = vp->z + localp->org.x;
  375.       globalp->dir.x = -localp->dir.y;
  376.       globalp->dir.y = -localp->dir.z;
  377.       globalp->dir.z = localp->dir.x;
  378.       break;
  379.     case WEST_WALL:
  380.       globalp->org.x = vp->x + localp->org.z;
  381.       globalp->org.y = vp->y - localp->org.y;
  382.       globalp->org.z = vp->z + localp->org.x;
  383.       globalp->dir.x = localp->dir.z;
  384.       globalp->dir.y = -localp->dir.y;
  385.       globalp->dir.z = localp->dir.x;
  386.       break;
  387.     default:
  388.       globalp->org.x = vp->x + localp->org.y;
  389.       globalp->org.y = vp->y + localp->org.x;
  390.       globalp->org.z = vp->z - localp->org.z;
  391.       globalp->dir.x = localp->dir.y;
  392.       globalp->dir.y = localp->dir.x;
  393.       globalp->dir.z = -localp->dir.z;
  394.       break;
  395.   }
  396. }
  397.  
  398. static void display_surface(void)       /* Display surfaces     */
  399. {
  400.   int col;              /* Element column index                 */
  401.   int row;              /* Element row index                    */
  402.   int surf;             /* Surface counter                      */
  403.   ELEM *elemp;          /* Element pointer                      */
  404.   SURF *surfp;          /* Surface pointer                      */
  405.   static int loop = 0;  /* Iteration counter                    */
  406.  
  407.   printf("Loop %d\n", loop++);
  408.   for (surf = 0; surf < NUM_SURF; surf++) {
  409.     printf("Surface %d\n", surf);
  410.     surfp = &(RoomSurf[surf]);
  411.     for (row = 0; row < surfp->num_row; row++) {
  412.       for (col = 0; col < surfp->num_col; col++) {
  413.         elemp = surfp->elemp + row * surfp->num_col + col;
  414.         printf("% 8.2f ", elemp->total / (PI * surfp->row_dim *
  415.             surfp->col_dim));
  416.       }
  417.       putchar('\n');
  418.     }
  419.   }
  420. }
  421.  
  422.